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ABSTRACT 



Aims. We present a general method to solve radiative transfer problems including scattering in the continuum as well as in lines in 
3D configurations with periodic boundary conditions. 

Methods. The scattering problem for line transfer is solved via means of an operator splitting (OS) technique. The formal solution is 
based on a full characteristics method. The approximate A operator is constructed considering nearest neighbors exactly. The code is 
parallelized over both wavelength and solid angle using the MPI library. 

Results. We present the results of several test cases with different values of the thermalization parameter and two choices for the 
temperature structure. The results are directly compared to ID plane parallel tests. The 3D results agree very well with the well-tested 
ID calculations. 

Conclusions. Advances in modern computers will make realistic 3D radiative transfer calculations possible in the near future. Our 
current code scales to very large numbers of processors, but requires larger memory per processor at high spatial resolution. 
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1. Introduction 



2. Method 
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Interest in 3-D radiative transfer in stellar atmospheres has 
grown with the calculations of Asplund and collaborators 



dAsplund et aljll999ll2000tlAsplundlEoorjl:lAsplund et a"Dl2005 



iGrevesse et al.l2007l) . This work has indicated that the solar oxy- 
gen abundance needs to be revised downward. However, the 
revised abundances are diffi cult to reco ncile with helioseismo- 
logical results (see Bas u & Antial 20081 and references therein). 
The work of Asplund et al. is based on comparisons of syn- 
thetic spectra produced by formal solutions of hydrodynamical 
models of solar convection. We present a framework for solv- 
ing the full scattering problem that is ap plicable to hydrody- 
namic al calculations of stellar atmospheres. [H auschildt & Baron 
(2006, hereafter: Paper I) and lBaron & Hausch ildt (2007[, here- 
after: Paper II) described a framework for the solution of the 
radiative transfer equation for scattering continua and lines in 
3D (when we say 3D we mean three spatial dimensions, plus 
three momentum dimensions) for the time independent, static 
case. In the 3rd paper of this series we apply these methods 
to problems with period boundary conditions which typically 
arise in radiation-hydrodynamical simulations of convective at- 
mospheres. In such calculations the radiation transport has to 
be simplified compared to the full problem in order to keep the 
calculations tractable. However, a full solution of the scattering 
line problem is needed for comparison and post-processing of 
the structures. 

We describe our method, its rate of convergence, and present 
comparisons to our well-tested 1-D calculations. 



In the following discussion we use notation of Papers I and II. 
The basic framework and the methods used for the formal so- 
lution and the solution of the scattering problem via operator 
splitting are discussed in detail in Papers I and II and will thus 
not be repeated here. 

In the following we assume (without restriction) that we 
have periodic boundary conditions in the x and y coordinates, 
and for the z coordinate that the 'bottom' (large optical depth) 
is at z = Zmin and the 'top' (interface to empty space) is at 
Z = Zmax- The implementation of the periodic boundary condi- 
tions within our framework is simple: We use a "full character- 
istics" approach that completely tracks a set of characteristics of 
the radiative transfer equation from the outer boundary through 
the computational domain to their exit voxel and takes care that 
each voxel is hit by at least one characteristic per solid angle. 
One characteristic is started on each boundary voxel (in this case 
these are the planes z = z m in and z = Zmax) and then tracked un- 
til it leaves the other boundary. The direction of a bundle of full 
characteristics is determined by a set of solid angles (6, <p) which 
correspond to a normalized momentum space vector (p x , p y , p z ). 
The periodic boundary conditions are simply implemented as 
a wrap-around (e.g., passing x max for p x > wraps around to 
x m i n ) and continuing of the characteristic until it leaves at the z 
boundary. Characteristics with very small \p z \ would require a 
large number of wrap-arounds (and eventually would lead to in- 
finitely long characteristics), therefore, we limit the number of 
wrap-arounds per voxel to a prescribed value, typically around 
16 (tests have shown that larger values do not affect the results, 
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values as small as 4 are usable in plane-parallel tests). The code 
is parallelized as described in Paper II. 

3. Plane-parallel tests 

3.1. Testing environment 

We use the framework discussed in Paper I and II as the baseline 
for the line transfer problems discussed in this paper. Our basic 
setup is similar to that discussed in Paper II. Periodic boundary 
conditions (PBCs) are realized in a plane parallel slab. We use 
PBCs on the x and y axes, z max is at the outside boundary, z m in 
the inside boundary. The slab has a finite optical depth in the z 
axis. The grey continuum opacity is parameterized by a power 
law in the continuum optical depth T st( j in the z axis. The basic 
model parameters are 

1. thickness of the slab, z max - z m in = 10 7 cm 

2. Minimum optical depth in the continuum, T™j n = 10~ 8 and 
maximum optical depth in the continuum, r™j X = 10. 

3. Constant temperatures (in all axes), T = 10 4 K 

4. Outer boundary condition, /r = and diffusion inner bound- 
ary condition for all wavelengths. 

5. Parameterized coherent & isotropic continuum scattering by 
defining 

Xc = e c K c + (1 - e c )cr c 

with < €c < 1 . k c and cr c are the continuum absorption and 
scattering coefficients. 

The line of the simple 2-level model atom is parameterized 
by the ratio of the profile averaged line opacity xi to the con- 
tinuum opacity Xc an d the line thermalization parameter e/. For 
the test cases presented below, we have used e c = 1 and a con- 
stant temperature and thus a constant thermal part of the source 
function for simplicity (and to save computing time) and set 
XilXc - 10 6 to simulate a strong line, with varying e\ (see be- 
low). With this setup, the optical depths as seen in the line range 
from 1(T 2 to 10 6 . We use 32 wavelength points to model the full 
line profile, including wavelengths outside the line for the con- 
tinuum. We did not require the line to thermalize at the center of 
the test configurations, this is a typical situation one encounters 
in a full 3D configurations as the location (or even existence) of 
the thermalization depths becomes more ambiguous than in the 
ID case. 

The slab is mapped onto a Cartesian grid. For the test cal- 
culations we use voxel grids with the same number of spatial 
points in each direction (see below). The solid angle space was 
discretized in (6, 0) with rig = if not stated otherwise. In the 
following we discuss the results of various tests. In all tests we 
use the full characteristic method for the 3D RT solution as de- 
scribed above. Unless otherwise stated, the tests were run on 
parallel computers using 128 CPUs. For the 3D solver we use 
n x = riy = n z = 2 * 32 + 1 = 65 points along each axis. The solid 
angle space discretization uses rig = = 64 points. 

3.2. Results 

We test the accuracy of the 3D PBC solution by comparing it 
to the results of the ID code for several line scattering param- 



eters. The ID solver uses 64 depth points, distributed logarith- 
mically in optical depth. Figures [T]43 show the mean intensities 
/ at T st d = and the z component of the emergent flux F as 
function of wavelength for both the ID (+ symbols) and the 3D 
solver. The agreement is excellent for all values of e/ from unity 
to 10~ 8 , indicating that the 3D code produces an accurate solu- 
tion even for extreme cases of line scattering. In the case with 
e/ = 1CT 8 the continuum processes lead to earler thermalization 
than the classical approximation J cc e 1 / 2 as the line strength is 
limited compared to the continuum. This behavior is the same as 
in the ID plane-parallel comparison case. The convergence rate 
of the line source function (here used together with Ng accelera- 
tion) is the same as discussed in Paper II, in the case of e/ = 1CT 8 
the 3D code needed 29 iterations with the nearest-neighbor A* to 
reach a relative accuracy of 10~ 8 using the simple starting guess 
S = B. The nearest-neighbor A* does allow stopping the it- 
erations earlier than a diagonal (local) A* due to the improved 
convergence rate (see paper I). This can easily cut the number 
of iterations by factors of two or more, even greater savings are 
possible if the accuracy limit is relaxed. 

In addition to the mean intensities, we checked that the flux 
vectors F have vanishing components in the x and y directions, 
typically maxflFJ, |F y |)/|F z | < 10~ 13 in all voxels. We stress that 
this result is the result of the calculations and is not forced by 
the numerical scheme. 



4. Tests with 3D structures 

For a test with a computed 3D structure, we obtained an exam- 



ple snapshot structure from H-G. Ludwig dCaffau et al.l 12007 



Wedemeyer et al. 2004) of a radiation-hydrodynamical Simula 
tion of convection in the solar atmosphere. The radiation trans- 
port calculations were performed with a total ofl41xl41xl51 
grid points in x,y, and z, respectively, for a total of 3 002 03 1 
voxels, the periodic boundary conditions were set in the (hori- 
zontal) x,y plane. The transport equation is solved for rig = 16 
and = 32 solid angle points, so that a total of about 1.5 x 10 9 
intensities are calculated for each iteration and wavelength point. 
For the tests described here, we are only using the temperature- 
density structure of the hydro model and ignore the velocity field 
for the simple tests presented here. We set the continuum opac- 
ity proportional to the density p by choosing a rough tempera- 
ture independent estimate for the Rosseland mean opacity per 
unit mass of 0.1 cm 2 /g and parameterizing the line in the same 
form as discussed above and in Paper II with a line total (wave- 
length integrated) opacity of 100 times the local continuum and 
32 wavelength points distributed over the full line profile. 

4.1. Results 

We ran a number of line transfer tests with ei = 1, 10~ 2 and 10~ 4 . 
The convergence rates for the two scattering cases are shown in 
Fig. [5] together with the convergence rates for a small plane par- 
allel test model and the results of the A iteration for continuum 
transfer (for computer time reasons) for a continuum e c = 10 2 . 
The convergence rate for the plane-parallel tests and the hydro 
model are remarkably similar, the A* operator delivers very rea- 
sonable and practically usable convergence rates. 
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4.1.1. Images 

Figures[6]and[7]show visualizations of the results for 3D contin- 
uum transfer. The RT problem was solved for e = 1 (left panels) 
and 10~ 2 (right panels) and a formal solution with the converged 
source functions was computed for given viewing angles. The 
graphs are actual images of the intensities as they would be seen 
by an external observer different angles. The visible surface is to 
the left, the 'sides' of the computational box could not be seen 
by an observer and are shown for information only. The effect 
of scattering on the images is similar to terrestrial fog in that it 
reduces the contrast of visible features; even moderate scattering 
of e c = 10 2 significantly reduces visibility. The limb darkening 
is also clearly visible in the figures. 

Figures [9] to [12] show images generated for the results of the 
line transfer solution. Three panels show results for individual 
wavelength (continuum, line wing and line center) and a com- 
posite image. The images are significantly different for these 
wavelengths. The line scattering produces a similar 'fog effect' 
as the scattering in the continuum transfer model, however, the 
images appear not that different. While one might expect that the 
line images would look vastly different from the continuum vi- 
sualization, part of the similarity is due to the fact that they were 
scaled individually in order to highlight the differences in struc- 
ture between the wavelengths rather than comparing them on a 
absolute scale. The composite image (best viewed in color avail- 
able in the online version of this paper) shows the differences in 
the visible structures between the different wavelengths. 

4.1.2. Limb darkening and contrast 

In Fig. [8] we show the limb darkening and contrast for a contin- 
uum test case with different values of e. To compute the limb 
darkening, we calculate the intensity average < / > over the vis- 
ible surface for different values of cos(6>) where 9 is the angle 
between the observer and the normal to the surface. We sim- 
ilarly calculate the contrast as (/- < / >) 2 >)/ < I > over 
the visible surface for different 9. The absolute values of the limb 
darkening and the contrast depend strongly on e, scattering dra- 
matically reduces the contrast and 'flattens' the limb darkening 
law. Overall, the limb darkening is nearly linear in cos(6>), as 
would be expected from a plane parallel atmosphere with grey 
temperature structure. 
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5. Conclusions 

Using rather difficult plane parallel test problems, we have 
shown that our 3D full-characteristics method gives very good 
results when compared to our well-tested ID code. The periodic 
boundary conditions method discussed here is particularly well- 
suited to 3-D hydrodynamical simulations of convection in stel- 
lar atmospheres and in future work we will compare our results 
to observations as well as to previous calculations. The results 
for the computed 3D structure show that A* leads also to good 
convergence for a true 3D structure, with convergence rates that 
are comparable to the simple test cases (see also Papers I and II). 

Acknowledgements. This work was supported in part by by NASA grants 
NAG5-3505 and NAG5-12127, NSF grants AST-0307323, and AST-0707704, 
and US DOE Grant DE-FG02-07ER41517, as well as DFG GrK 1351 and 



Hauschildt and Baron: 3D radiative transfer framework III 



10' 1 
9.900X10 3 



ine transfer, epsilon(line) = 1 e — 



9.950X10 3 1.000x10* 1.005x10* 1.010x10* 

wavelength 

ine transfer, epsilon(line) = 1e— 



ne transfer, epsilon(line) = 1 e — 4 




E 9.900X10 3 9.950X10 3 1.000x10* 1.005x10* 1.010x10* 

wavelength 

_ line transfer, epsilon(line) = 1e — 4 



10" 
9.900x1 3 



1.000x10* 1.005x10 
wavelength 



10" 




9.900x10 3 9.950X10 3 



1.000x10 
wavelength 



Fig. 1. The mean intensity J and the z component of the radia- 
tion flux F at r st( j = as function of wavelength. The + symbols 
are the comparison results with the ID solver, the full lines the 
results from the 3D PBC solution. The results are for e\ — 1 and 
constant temperatures. 
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Fig. 2. The mean intensity J and the z component of the radia- 
tion flux F at r st( j = as function of wavelength. The + symbols 
are the comparison results with the ID solver, the full lines the 
results from the 3D PBC solution. The results are for e; = 10~ 2 
and constant temperatures. 



Fig. 3. The mean intensity J and the z component of the radia- 
tion flux F at Tjta = as function of wavelength. The + symbols 
are the comparison results with the ID solver, the full lines the 
results from the 3D PBC solution. The results are for e/ = 10 -4 
and constant temperatures. 
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Fig. 4. The mean intensity J and the z component of the radia- 
tion flux F at Tjtd = as function of wavelength. The + symbols 
are the comparison results with the ID solver, the full lines the 
results from the 3D PBC solution. The results are for e/ = 10 
and constant temperatures. 
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Fig. 8. Limb darkening (left panel) and contrast ^((1 - (/)) 2 ))/ (/) for 3D continuum transfer with the hydro structure. 
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Fig. 5. Convergence rates of the 3D transfer for line transfer 
with plane-parallel test structures (label 'PP') and the 3D hydro 
structure (label 'hydro'). For comparison, the convergence of the 
A iteration for plane-parallel continuum transfer is also shown. 



Fig. 9. Visualization of the results for the line 3D radiation trans- 
fer with e; = 1. The images are intensities in the directions 
<p = 25 deg and 8 = deg. The top left panel is the image in 
the continuum, the top right panel the image at the line center, 
the bottom left panel the image in the line wing, the bottom right 
panel is a composite image. 



Fig. 10. Visualization of the results for the line 3D radiation 
transfer with e/ = 1CT 4 . The images are intensities in the direc- 
tions cp — 25 deg and 8 = deg. The top left panel is the image 
in the continuum, the top right panel the image at the line center, 
the bottom left panel the image in the line wing, the bottom right 
panel is a composite image. 



Fig. 11. Visualization of the results for the line 3D radiation 
transfer with e\ — 1 . The images are intensities in the directions 
(f> = 25 deg and = 50 deg. The top left panel is the image in 
the continuum, the top right panel the image at the line center, 
the bottom left panel the image in the line wing, the bottom right 
panel is a composite image. 



Fig. 6. Visualization of the results for continuum 3D radiation 
transfer for e c — 1 (left panel) and 10~ 2 (right panel). The images 
are intensities in the directions <p — 25 deg and 9 = deg (top 
row) and 8 = 40 deg (bottom row). An observer would only see 
the left face of the cube (inside the indicated area), the other sides 
of the cube are shown for clarity and are actually invisible due to 
the periodic boundary conditions. The scaling of the intensities 
is the same within each column but different for the left and right 
columns. 

Fig. 7. Same as Fig. [6] but for 8 = 60 deg (top row) and 8 = 
80 deg (bottom row). The scaling of the intensities is as in Fig. [6] 
The effect of limb darkening is clearly visible in this figure. 



Fig. 12. Visualization of the results for the line 3D radiation 
transfer with e/ = 10 -4 . The images are intensities in the direc- 
tions <p — 25 deg and 8 = 50 deg. The top left panel is the image 
in the continuum, the top right panel the image at the line center, 
the bottom left panel the image in the line wing, the bottom right 
panel is a composite image. 
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